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Abstract 

Restricting the randomization of hard-to-change factors in industrial experiments is often 
performed by employing a split-plot design structure. From an economic perspective, these 
designs minimize the experimental cost by reducing the number of resets of the hard-to- 
change factors. In this paper, unbalanced designs are considered for cases where the subplots 
are relatively expensive and the experimental apparatus accommodates an unequal number 
of runs per whole-plot. We provide construction methods for unbalanced second-order split- 
plot designs that possess the equivalence estimation optimality property, providing best linear 
unbiased estimates of the parameters; independent of the variance components. Unbalanced 
versions of the central composite and Box-Behnken designs are developed. For cases where 
the subplot cost approaches the whole-plot cost, minimal point designs are proposed and 
illustrated with a split-plot Notz design. 
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Introduction 


Frequently, industrial experiments involve factors that are difficult, time consuming, or costly 
to manipulate, referred to as hard-to-change (HTC) factors. Alternatively, some factors may 
be relatively easy to manipulate, referred to as easy-to- change (ETC) factors. A split-plot 
design structure is an efficient experimental approach that reduces the number of settings of 
the HTC factors. We assign the HTC factors to the whole-plots and the ETC factors to the 
subplots. For a more detailed development of split-plot designs see Montgomery (2004) and 
Myers and Montgomery (2002). 

Vining, Kowalski, and Montgomery (2005) and Parker, Kowalski, and Vining (2005a) pro- 
posed design construction techniques for balanced equivalent estimation designs. To achieve 
balance (same number of subplot runs in each whole-plot) in these designs, subplot runs 
were replicated within each whole-plot. Unbalanced versions, with different whole-plot sizes, 
reduce the replication of the subplot runs allowing a more efficient allocation of experimental 
resources at the subplot level. In addition, when the whole-plot size is large, replicating sub- 
plot runs to achieve balance generates an excessive number of subplot error variance degrees 
of freedom. Unbalanced designs provide attractive alternatives to practitioners when they 
are permissible. 

In experimental design, we strive for small designs that meet our objectives and simulta- 
neously minimize cost. Similarly in the split-plot context, the cost of the experiment is 
an important criteria, however the definition of small is not as clear because we have two 
types of experimental units with different associated costs. In many cases, the cost of the 
whole-plot experimental units is much greater than that of the subplot units, making the 
subplot runs essentially free. A balanced design is appropriate in this situation. However, if 
the subplot runs are relatively expensive as compared to subplot runs, then an unbalanced 
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design that simultaneously minimizes the number of whole-plots and subplot experimental 
units is attractive from an economic perspective. Note that if the cost of the whole-plots 
equals the cost of the subplots, then a completely randomized design is appropriate. 

In addition to the cost considerations, the experimental apparatus plays a role in determining 
whether an unbalanced design is a viable option. For example, if the subplot units are 
specimens and the whole-plots represent the settings of a temperature chamber, then varying 
the number of subplot experimental units may introduce an additional source of variability 
due to changing the thermal load in the chamber. 

In aerodynamic wind tunnel research, the configuration of the scaled test article usually re- 
quires a partial mechanical disassembly, making it a HTC factor. Alternatively, the flow held 
parameters are relatively ETC change (see Parker, Kowalski, and Vining (2005a, 2005b)). 
To simulate flight conditions, there is a class of wind tunnel laboratories that operate in 
a cryogenic pressurized environment. In these complex facilities, the setting of a whole- 
plot factor condition is time consuming and expensive. This is due to the conditioning of 
the environment that allows personnel access the test article and perform the configuration 
change. Once the whole-plot conditions are set, maintaining these extreme environmental 
conditions over time is costly and is directly proportional to the number of subplot runs. In 
this example, the whole-plot size is not determined by the capacity of the experimental ap- 
paratus, rather it represents the time to execute the subplot runs. An economical design for 
this application simultaneously minimizes the number of whole-plots and the total number 
of subplot runs. 

In this paper, we review the conditions under which ordinary least squares is an appropriate 
method of estimation for a second-order split-plot design. Based on these conditions, we 
illustrate two systematic design construction strategies to build unbalanced equivalent esti- 
mation designs from the central composite (CCD) and Box-Bchnken (BBD). Replication of 
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subplot runs within each whole-plot has been reduced over the balanced versions, thereby 
reducing the size of some whole-plots. This development does not attempt to generalize all 
of the admissible forms of unbalance. Instead, we focus on a strategy to reduce the size of 
the whole-plots containing subplot replicates, especially those containing replicated subplot 
centers. In addition, we discuss minimizing the total number of subplot runs by configuring 
response surface minimal point designs into an equivalent estimation split-plot structure. 

The notation D (w,k) is used to denote a design with w whole-plot (HTC) factors and k 
subplot (ETC) factors. We use a and f3 to denote the distance to the axial points for the 
subplot and whole-plot factors, respectively. We assume that the value of a is the same for 
all of the subplot factors and similarly /3 is the same for all of the whole-plot factors. Designs 
with combinations of 1 < w < 3 and 2 < k < 4 are considered to encompass many practical 
situations. Note that unbalanced designs with k = 1 result in whole-plots of size one and 
have not been included. 

Second-Order Split-Plot Designs 

Letsinger, Myers, and Lentner (1996) emphasized the need for research in the area of second- 
order split-plot designs and focused on analysis issues; recommending the use of restricted 
maximum likelihood for model estimation and inference. They proposed sorting a completely 
randomized design (CRD) into a split-plot structure, which often results in an unbalanced de- 
sign. Draper and John (1998) considered modification of the central composite design (CCD) 
(Box and Wilson (1951)) to be executed in a split-plot mode to achieve near-rotatability. 
Trinca and Gilmour (2001) studied multiple layers of split-plotting, called multistratum de- 
signs, where the two-level split-plot design is a special case of this general class. They 
employed a computer-intensive search algorithm to construct designs that maintained near- 
orthogonality between strata. Goos and Vandebroek (2001, 2003) proposed the use of a 
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point-exchange algorithm to construct ©-optimal split-plot designs, which minimize the vol- 
ume of the joint confidence region of regression coefficients. Goos (2002) provides a thorough 
development of balanced and unbalanced ©-optimal split-plot designs. 

Vining, Kowalski, and Montgomery (2005), hereafter referred to as VKM, proposed a class 
of second-order SPD’s where ordinary least squares (OLS) is equivalent to generalized least 
squares (GLS) for model estimation. Moreover, they proposed augmentation of these designs 
to provide pure-error estimates at both the whole-plot and subplot levels. Parker, Kowalski, 
and Vining (2005a) provided a generalized derivation of the equivalence conditions that leads 
to the development of systematic construction strategies and the ability to numerically verify 
if an arbitrary design achieves the equivalence property. They illustrated the construction 
of split-plot versions of the central composite design and a class of equivalent estimation 
©-optimal designs. Parker, Kowalski, and Vining (2005b) extends the classes of designs to 
include balanced split-plot versions of the Box-Behnken, equiradial, and small composite 
design. 

Equivalent estimation designs offer a number of attractive features. For example, the equiv- 
alence property is independent of the variance components. Moreover, the parameter esti- 
mates are BLUE (best linear unbiased estimators), independent of the variance components, 
and robust to the assumption of normality. Model estimation of equivalent estimation de- 
signs is simplified using OLS; eliminating the need for specialized software that is often times 
unavailable to the industrial practitioner. 

Equivalence Conditions 

We define the general form of the model as 

y = X/3 + S + e , 
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where y is the N x 1 vector of responses, X is the N x p model matrix with rank of p, (3 is 
the p x 1 vector of coefficients, S is the N x 1 vector of random whole-plot errors, e is the 
JV x 1 vector of random subplot errors, N is the total number of subplot runs, and p is the 
number of terms in the model including the intercept. We assume that S + e has a mean 
of 0 and variance-covariance matrix X = erf I + crj J&, where of is the subplot error variance, 
and erf is the whole-plot error variance. For an unbalanced design, J& is a function of the 
N x p incidence matrix Z expressed as J& = ZZ' 

The matrix J& has the form of 

'Win, 0 0 

T o i n2 i; 2 ••• o 

j b — 

0 

where m is the number of whole-plots and n % is the number of subplot runs in the i th whole- 
plot and denotes the length of the vector of ones. 

Parker, Kowalski, and Vining (2005a) provide a detailed derivation of the equivalence con- 
dition, summarized as follows. We assume that there is one level of split-plotting and the 
design supports the specified model. The necessary and sufficient condition for equivalence 
from McElroy (1967) (see also Graybill (1976, p. 209)) is 

XF = XX, 

where F is a p x p non-singular matrix. We define F to be 

F = erf I + crfK, 


1 1 ' 

J -n r n- L n rl 


where, 

K = (X'X) -1 X'J fe X. (1) 
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It can be shown that the condition for equivalence can be expressed as 


XK = J 6 X 


( 2 ) 


The right-hand side of this equation sums the columns of the model matrix within each 
whole-plot. These column totals, or unsealed whole-plot means, must equal the model matrix 
multiplied by K to satisfy the condition of equivalence. Note that the equivalence condition 
is independent of the variance components, erf and erf. 


For a complete second-order model, let X be partitioned as 



1 

w Dl 

W Ql 

s Dl 

s Ql 

X = 

1 

W Da 

Wq 2 

Sd 2 

Sq 2 


1 

w Dm 

w Qm 

s Dm 

Sq 


where the first column corresponds to the intercept, Wd ; denotes the whole-plot main ef- 
fects and two-factor interactions, Wq contains the whole-plot pure quadratics, S D; includes 
the subplot main effects, subplot by subplot two-factor interactions, and the whole-plot by 
subplot interactions, and Sq. denotes the subplot pure quadratic terms of the i th whole-plot. 


Consider a form of K for the unbalanced design as, 


n 0 

O' 

0' 

O' 

m n 

0 

uwId 

0 

0 

0 

u w 

0 

nwlq 

0 

Mw 

0 

0 

0 

0 

0 

u s 

0 

0 

0 

Vg 


( 3 ) 


where n 0 and % are scalars, m 0 is a vector of length k, Id and Iq are (w + U) C’- 1 ) ) x (w + 
) and w x w identity matrices, respectively, uw is a vector of length w, us is a vector 
of length k, Mw is a w x k matrix, and Vs is a k x k matrix. The definition of m 0 , Mw, 
and Vg were defined by Parker, Kowalski, and Vining (2005a, 2005b) for balanced designs. 
The proposed elements, uw, and ug, are specific to unbalanced designs and are defined in 
this paper. 
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Performing the multiplication in Equation (2) we obtain 


1 w Dl w Ql s Dl s Ql 


XK = 


1 W, 


Q 2 °D2 


n 0 O' 

0 uwId 


O' 0' m' 0 

0 0 0 


0 uwIq 0 Mv 
0 0 0 0 


1 W D Wo 

s Dm s 

Qm J [ Ug 

0 

0 

nil 

ll'W Dl 

ll'W Ql 

ll'S Dl 

ll'Sq, 

n 2 l 

J b x = ; 

ll'W D2 

H'Wq 2 

ll'S D2 

H'Sq 2 

n m 1 

ll'W Dm 

II'Wq™ 

11'Sd™ 

11'Sq. 


where the vectors of ones have length n t in the i th whole-plot (explicit subscripts on these 
vectors are omitted for clarity). From Equation (4), we see that the conditions to achieve 
equivalence are isolated to the relationships between the within whole-plot column sums of 
the intercept and the whole-plot and subplot pure quadratic terms. In addition, we see by 
inspection that the designs possess first-order plus interaction orthogonal subplot designs. 

If the design possess the equivalence property, then /3 0 ls = @gls an d ordinary least squares 
(OLS) is an appropriate method of estimation. Note that it is a rather remarkable result that 
OLS would be an appropriate method of estimation for a class of unbalanced second-order 
split-plot designs. The OLS estimate of the model coefficients is f3 OLS = (X'X) _ 1 X'y. This 
expression highlights that the parameter estimates from an equivalent estimation split-plot 
design are independent of the variance components. By equivalence to GLS, the variance- 
covariance matrix for the OLS estimate is (X'E _ 1 X) _1 . 


Balanced Design Summary 


Parker, Kowalski, and Vining (2005a) proposed two systematic design construction tech- 
niques that enable the construction of equivalent estimation SPD’s derived from completely 
randomized response surface designs. The first technique is a generalized version of the 



Table 1: Elements of K for Balanced Designs. 

(VKM and MWP denotes method of construction, Vs, mo, and 0 denote the elements of K) 



VKM-CCD 

VKM-BBD 

MWP-BBD 

V s 

n -I -I / 

2(*-i)i*n 

0 

m 0 

0 

0 

4(& — l)lfc 

4> 

0 

0 

6 — 4 k 


construction method proposed by Vining, Kowalski, and Montgomery (2005), subsequently 
referred to as the VKM method. In the second method, we obtain designs that achieve the 
minimum number of whole-plots to configure a classical CCD into a split-plot structure, 
subsequently referred to as the MWP method. Based on the construction method and the 
parent completely randomized design (CRD), we define the elements of K. 

The elements of K in Equation (3), derived in Parker, Kowalski, and Vining (2005a, 2005b) 
are summarized in Table 1. For the balanced case, Uq = % = n and is the number of 
subplot runs per whole-plot. We define Mw to be of the form 01„,l / fc . Note that the values 
for the VKM construction method apply to cases where the number of subplot factors is a 
power of 2, and the values for the MWP method apply to w — 1. 

Unbalanced VKM Method CCD 

To construct an unbalanced design, we first construct the balanced version. Then, we reduce 
the size of the whole-plot containing all factors at their center level. We denote the size of 
this whole-plot as n c in the derivation, and choose n c = 2 in the examples presented. 

The construction of a VKM method CCD is described in Parker, Kowalski, and Vining 
(2005a). We begin with a completely randomized CCD in w + k — f factors. The design 
is rearranged such that the HTC factors remain constant within each WP and the subplot 
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designs are orthogonal for a first-order pins interaction model. The axial points (or star 
points) of the HTC factors are placed in a WP by themselves. The center points, consisting 
of all factors at their zero level in coded units, are also placed in a separate whole-plot. To 
maintain balance, subplot runs are replicated as required within each whole-plot. 

For an unbalanced VKM method CCD, when k is a power of 2, we define K by setting 
mo = 0 and Mw = 0. Based on this form of K and Equation 4, the conditions that must 
be satisfied are 


n w W Di - l ni l( l .W Di 

»wW Qj = DI^Wq. V i (5) 

and 

W'A = S Qi V s 

riilm = n 0 l ni + W Q .u w + S Q .u s Vi. (6) 

For the conditions in Equation (5), since whole-plot factors are held constant within each 
whole-plot, the column sums of the whole-plot model matrix equals n t times the level of 
the whole-plot model term. Therefore, we End % equals the number of subplot runs in 
the whole-plots containing the HTC factors at their non-zero levels. Note that for this 
construction method all whole-plots of this type must possess the same number of subplot 
runs. These whole-plot model conditions are satisfied for any design by the nature of the 
split-plot structure. 

In Equation (6), the first condition involves the column sums of the subplot pure quadratic 
terms and Vg is defined in Table 1. The second condition expresses a relationship between 
each row of the whole-plot model matrix and the number of subplot runs in that whole-plot, 
n t . We denote the j th row of the model matrix for the whole-plot and subplot pure quadratic 
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terms within the i th whole-plot as Wq and Sq , respectively. We define the general form 
of u w to be Al^ and u s to be pl k , where A and p are suitable constants chosen to satisfy 
the equivalence conditions. Therefore, we express the unbalanced condition as 

rii = n 0 + XW Qij l w + pS Q ..l k V i,j. (7) 

To find expressions for no, A, and p, we consider the forms of Wq and Sq. . for each type 
of whole-plot found in a CCD constructed according to the VKM method, summarized in 
Table 2. Note that for the whole-plot and subplot axial runs the position of a 2 and f3 2 will 
change depending on the factor involved, however there is always a single element in each 
row vector. 


Table 2: VKM Method CCD Cases of Whole-Plot Type. 

(Wq,, is the j th row of the model matrix containing the WP pure quadratics, 
Sq is the j th row of the model matrix containing the SP pure quadratics, 
n t is the number of subplot runs in the i th whole-plot, 
a and f3 are the axial point distances for the ETC and HTC factors, respectively) 


Case 

WP Type 

Wq 17 

s q. 

n t 

c 

all centers 

[o 

■■■ °L 

i o ■■■ 

°L 

n c 

1 

WP centers, SP axials 

[0 

••• 0] 

J w 

[ a 2 0 ■ ■ 

• 0]* 

rii = 2A: 

2 

WP axials, SP centers 

[ip 

0 ••• 0 1 

J w 

[o 


n 2 — n f 

3 

factorial 

[1 

... 1] 
ilil 

ri - 


CO 

S 


From case (c) and Equation (7), we see that no = n c . Based on the structure of the row 
vectors, we can express Equation (7) for each case as 

Case (1): n\ — n c = pa 2 

Case (2): n 2 — n c = X/3 2 
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Case (3): 713 — n, 


A w + pk. 


In a CCD, ri\ = 2k and ri2 = 713 = rif, where rif denotes the number of subplot runs in the 
factorial whole-plots. Therefore, we find 

2k — n c 


n f — n c 
x = ^- 


In addition from case (2) and (3), we see that A/3 2 = A w + pk. Solving this expression for a 
and (3 provides, 


a = 


k(2k — n c 


P = 


^ {nf - n c ) [1-^ 
wa 2 (nf — n c ) 


a 2 [rif — n c ) — k(2k — n c ) 
If we desire a design with a = /3, then we obtain 


a = P = 4 w + 


k(2k — n c ) 
i n f - n c ) 


( 8 ) 

(9) 


Performing an analysis of Equations (8) and (9), we find that our choice of a and P are 
constrained by 


a > 


'k{2k — n c 
rif - n c 


P > \/w . 


In addition from Equations (8) and (9), we see that the relationship between a and P is 
inversely proportional. Using these expressions, we can define K for an unbalanced VKM 
CCD when k is a power of 2. 
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Table 3: VKM Method Unbalanced Central Composite Design for D(2,2). 
and z 2 are HTC factors, Xi and x 2 are ETC factors, a and (3 are the distances 
to the axial points for the ETC and HTC factors, respectively) 


Whole-Plot 

Zl 

Z2 

X\ 

X2 

Whole-Plot 

Zl 

Z"2 

X\ 


1 

-1 

-1 

-1 

-1 

6 

P 

0 

0 

0 


-1 

-1 

1 

-1 


(3 

0 

0 

0 


-1 

-1 

-1 

1 


(3 

0 

0 

0 


-1 

-1 

1 

1 


(3 
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0 

0 

2 

1 

-1 

-1 

-1 

7 

0 

-p 

0 

0 


1 

-1 

1 

-1 


0 

-(3 

0 

0 


1 

-1 

-1 

1 


0 

-(3 

0 

0 


1 

-1 

1 

1 


0 

-(3 

0 

0 

3 

-1 

1 

-1 

-1 

8 

0 

(3 

0 

0 


-1 

1 

1 

-1 


0 

(3 

0 

0 


-1 

1 

-1 

1 


0 

(3 

0 

0 


-1 

1 

1 

1 


0 

(3 

0 

0 

4 

1 

1 

-1 

-1 

9 

0 

0 

-a 

0 


1 

1 

1 

-1 


0 

0 

a 

0 


1 

1 

-1 

1 


0 

0 

0 

-a 


1 

1 

1 

1 


0 

0 

0 

a 

5 

-P 

0 

0 

0 

10 

0 

0 

0 

0 


-P 

0 

0 

0 


0 

0 

0 

0 


-P 

0 

0 

0 







-P 

0 

0 

0 







Consider a split-plot design with two hard-to-change factors, denoted by z i and z 2 , and 
two easy-to-change factors, denoted by X\ and x 2 . An unbalanced central composite design 
(CCD) in coded units is provided in Table 3. For this design, we have w — k — 2, n c — 2, 
rif = 4, and choose a — (3 — 2, resulting in n 0 = 2, % = 4, u w = \l W) and u s = \l k) 
and V s = flfcl'k = 21 fc l' fe , where n refers to the number of subplot runs in the whole-plot 
containing the subplot axial points. 

A balanced version of this design is identical except that whole-plot number 10 contains 
4 subplot runs. At first glance, the small reduction in design size may not appear to be 
significant. However, we suggest augmenting the base design by replication of whole-plot 


13 



number 10 to enable a pure-error based estimate of whole-plot error variance and to increase 
the degrees of freedom for a pure-error estimate of subplot variance. For df w degrees of 
freedom for the whole-plot error estimate, df w + 1 replicates of whole-plot 10 are required. 
Therefore, we see that the reduction in design size is not merely 2 subplot runs, rather it 
is 2 (df w + 1). The smaller size of these replicated whole-plots should be more appealing 
to a practitioner who is contemplating the cost-benefit trade-off of performing replicated 
whole-plots. This is especially true in the context that we are considering where the subplot 
runs are costly. 

When k is not a power of 2, an unbalanced design can be generated if we eliminate one set of 
whole-plots containing the whole-plot factor at its axial level. The resulting design retains 
the capability to estimate a complete second-order model, however it does not contain all 
the design points found in a classical central composite design. An attractive feature of this 
structure is that the distance to the axial points is unrestricted. For example, consider a 
design for D(l, 3) in Table 4. The form of K for this type of design requires a non-zero M\y, 
which is not found in a classical VKM construction. A general development of this form 
of K is not provided, however for the design in Table 4, K is defined by n 0 = 2, nyy = 8, 
uw = — 6 , us = 41*;, Vs = 21fcl' fc , m 0 = 0 , and Mw = 21' fc . Note that this design has two 
whole-plots with eight subplot runs, one with six subplot runs, and one with two subplot 
runs. 

A summary of unbalanced VKM CCD’s is provided in Table 5. The allotment of subplot runs 
to the whole-plots is denoted by a whole-plot size followed by its frequency in parentheses. 
For example, the design in Table 4 is denoted by 8(2), 6(1), 2(1). The D(w,3) entries are 
based on removing one set of whole-plot axial points, as previously discussed, and therefore 
require two fewer whole-plots than other designs for the same w. 
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Table 4: VKM Method Unbalanced Central Composite Design for D(l,3). 
(zi is a HTC factor, x\, x-z, and X 3 are ETC factors, a is the distance 
to the axial points for the ETC factors.) 


Whole-Plot 

Zl 

Xi 

X2 

%3 

1 

-1 

1 

-1 

-1 


-1 

-1 

1 

-1 


-1 

-1 

-1 

1 


-1 

1 

1 

1 


-1 

-1 

-1 

-1 


-1 

1 

1 

-1 


-1 

1 

-1 

1 


-1 

-1 

1 

1 

2 

1 

1 

-1 

-1 


1 

-1 

1 

-1 


1 

-1 

-1 

1 


1 

1 

1 

1 


1 

-1 

-1 

-1 


1 

1 

1 

-1 


1 

1 

-1 

1 


1 

-1 

1 

1 

3 

0 

-a 

0 

0 


0 

a 

0 

0 


0 

0 

-a 

0 


0 

0 

a 

0 


0 

0 

0 

-a 


0 

0 

0 

a 

4 

0 

0 

0 

0 
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0 

0 

0 
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Table 5: Summary of Unbalanced Central Composite and Box-Behnken Equivalent Estima- 
tion Split-Plot Designs by Construction Method. 


(w is the number of HTC factors, k is the number of ETC factors, m is the number of 
whole-plots, n is the number of subplot runs per whole-plot denoted by the whole-plot size 
with its frequency in parentheses, and N is the total number of subplot runs.) 


Factors 


VKM-CCD 



VKM-BBD 



MWP-BBD 


w 

k 

m 

n(freq.) 

N 

m 

n(freq.) 

N 

m 

n(freq.) 

N 

1 

2 

6 

4(5), 2(1) 

22 

4 

4(3), 2(1) 

14 

3 

5(1), 4(2) 

13 

1 

3 

4 

8(2), 6(1), 2(1) 

24 

4 

12(3), 2(1) 

38 

3 

13(1), 6(2) 

25 

1 

4 

6 

8(5), 2(1) 

42 

4 

24(3), 2(1) 

74 

3 

25(1), 8(2) 

41 


2 

2 

10 

4(9), 2(1) 

38 

10 

4(9), 2(1) 

38 

2 

3 

8 

6(1), 4(6), 2(1) 

32 

10 

12(9), 2(1) 

110 

2 

4 

10 

8(9), 2(1) 

74 

10 

8(9), 2(1) 

74 

3 

2 

16 

4(15), 2(1) 

62 

20 

4(19), 2(1) 

78 

3 

3 

14 

6(1), 4(12), 2(1) 

56 

19 

4(18), 2(1) 

74 

3 

4 

16 

8 (15), 2(1) 

122 

21 

4(20), 2(1) 

82 
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Unbalanced VKM Method BBD 


The Box-Behnken design is an attractive alternative to the central composite in the com- 
pletely randomized context due to its competitive size and utilization of three-levels. In the 
split-plot context, the three-levels are particularly appealing for the HTC factors. In this 
section, we consider applying the VKM construction method to the BBD. 

As a result of our construction philosophy to utilize the minimum number of whole-plots, the 
number of subplot runs per whole-plot increases with the number of subplot factors. These 
large whole-plots may not be practical in certain applications, thereby requiring a design 
with additional whole-plots of a smaller size. Nevertheless, in this section we assume that 
these minimum whole-plot designs are permissible and consider methods to reduce the size 
of the whole-plots containing replicated subplot runs. 

To minimize the number of whole-plots, we construct a single whole-plot with the HTC 
factor set at its center level. The size of this whole-plot, for / < 5, is 

2 ) = 2k(k — 1). 

To achieve balance, the subplot designs are replicated as required in the other whole-plots. 
The general form of K and the conditions expressed in Equation (6) hold for the BBD case. 
Following a similar approach as the CCD, we summarize the forms of Wq and Sq (j found 
in each type of whole-plot in Table 6. For cases (l)-(3), the position of the ones will change 
depending on the factors involved, however the form is consistent. Note that case (3) is not 
present when w — 1. 

From case (c) and Equation (7), we see that n 0 = n c . From case (1), we find 

n\ - n c = 2 p 

rti — n c 2 k{k — 1) — n c 
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Table 6: VKM Method BBD Whole-Plot Types. 


(Wq is the j th row of the model matrix containing the WP pure quadratics, 
Sq is the j th row of the model matrix containing the SP pure quadratics, 
n t is the number of subplot runs in the i th whole-plot) 


Case 

Wq, 

Sq, 

Hi 

c 

o 

o 

O 

O 

n c 

1 

[0 ••• 0] 

L J w 

O 

O 

?F 

ni = 4 ( 2 ) = - !) 

2 

O 

O 

S 

o 

o 

ri2 

3 

[1 1 o ••• o 1 

L J w 

O 

o 

n 3 


For case (2), we find 

ri ‘2 - n 0 = A + p 
A = n 2 — no — p . 

Similarly for case (3), we find 

2 k(k — 1) — n c 
A= 2 • 

Note that if ri\ = ri 2 — n 3 — 2k(k — 1 ), then A = p. 

Consider a design for the D(l,4) case given in Table 7. To construct this design, we begin 
with a balanced BBD consisting of 4 whole-plots of size 24. In whole-plots 1 and 2, there are 
3 replicates of the base subplot design to balance the design. The vector notation denotes 
two subplot runs for each row in whole-plots 1 and 2. As with the CCD, we reduce the size 
of the whole-plot containing all factors at their center level. The number of subplot runs is 
reduced by 22(df w + l), which is very significant due to the large whole-plot size. Moreover, it 
is apparent that the likelihood of running whole-plot replicates is improved considerably. We 
note that the extreme unbalancing in this design may not be desirable in many experimental 
circumstances, however it does illustrate the achievable reduction in subplot replication. 
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Table 7 : VKM Method Unbalar 
(zi is a HTC factor, xi , x 


Whole-Plot 

Z \ 

Xi 

%2 

X3 

1 

-1 

±1 

0 

0 


-1 

0 

±1 

0 


-1 

0 

0 

±1 


-1 

0 

0 

0 


-1 

±1 

0 

0 


-1 

0 

±1 

0 


-1 

0 

0 

±1 


-1 

0 

0 

0 


-1 

±1 

0 

0 


-1 

0 

±1 

0 


-1 

0 

0 

±1 


-1 

0 

0 

0 

2 

+1 

±1 

0 

0 


+1 

0 

±1 

0 


+1 

0 

0 

±1 


+1 

0 

0 

0 


+1 

±1 

0 

0 


+1 

0 

±1 

0 


+1 

0 

0 

±1 


+1 

0 

0 

0 


+1 

±1 

0 

0 


1 0 

0 ±1 
0 0 


, x 3 and £4 are ETC factors.) 


£4 

Whole-Plot 

Zl 

Xi 

%2 

X3 

X4 

0 

3 

0 

0 

-1 

-1 

0 

0 


0 

0 

1 

-1 

0 

0 


0 

0 

-1 

1 

0 

±1 


0 

0 

1 

1 

0 

0 


0 

-1 

0 

0 

-1 

0 


0 

1 

0 

0 

-1 

0 


0 

-1 

0 

0 

1 

±1 


0 

1 

0 

0 

1 

0 


0 

0 

0 

-1 

-1 

0 


0 

0 

0 

1 

-1 

0 


0 

0 

0 

-1 

1 

±1 


0 

0 

0 

1 

1 

0 


0 

-1 

-1 

0 

0 

0 


0 

1 

-1 

0 

0 

0 


0 

-1 

1 

0 

0 

±1 


0 

1 

1 

0 

0 

0 


0 

0 

-1 

0 

-1 

0 


0 

0 

1 

0 

-1 

0 


0 

0 

-1 

0 

1 

±1 


0 

0 

1 

0 

1 

0 


0 

-1 

0 

-1 

0 

0 


0 

1 

0 

-1 

0 

0 


0 

-1 

0 

1 

0 

±1 


0 

1 

0 

1 

0 


0 0 0 0 
0 0 0 0 


For the design in Table 7, we have w — 1, k — 4, n c — 2, resulting in no = 2, % = 24, 
u w = 11 , u s = (ll)l fc , V s = 61 fc l' fc , m o = 0, and M w = 0. 

A potentially unattractive feature of the design in Table 7 is the replicated subplot designs in 
whole-plots 1 and 2. We can reduce the size of these whole-plots by adding a non-zero Mw 
matrix to K. A general development of this form of K is not provided, however to illustrate 
this approach consider the design in Table 8. The first two whole-plots now contain 8 subplot 
runs each with no replication. For this design. K is defined by no = 2, % = 8, uw = —5, 
u s = (11)1*, V s = 61 fcl'j., M w = -41' fc , and m 0 = 0. 

The unbalanced design in Table 8 has achieved a significant reduction with m = 4 and 
N = A2 compared to Table 7 with rn — 4 and N = 70, and the balanced version, which 
requires rn — 4 and N = 96. Even though the number of whole-plots has not been reduced, 
the reduction in design size and particularly the size of the whole-plots that are proposed to 
be replicated, makes this design an attractive alternative. 

A summary table of unbalanced VKM BBD’s is provided in Table 5. For D(l,/c), these 
designs require the same or fewer whole-plots than their CCD competitors, however the 
larger whole-plot size may not be permissible in some applications. In contrast, the D(3 ,k) 
designs require additional whole-plots over the unbalanced CCD’s. 

Unbalanced MWP Method BBD 

In this section, we propose a class of unbalanced equivalent estimation designs that requires a 
minimum number of whole-plots to configure a completely randomized BBD into a split-plot 
structure for the w = 1 case. We begin with a balanced design constructed using the MWP 
method. This method does not require the subplot centers to be contained in a whole-plot by 
themselves, rather they can be included in other whole-plots and added one-by-one without 
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Table 8: Minimum Size VKM Method Unbalanced Box-Behnken Design for D(l,4). 
is a HTC factor, xi, x^ r x% and X 4 are ETC factors.) 


Whole-Plot 

Zl 

Xi 

X2 

X3 

X4 

1 

-1 

±1 

0 

0 

0 


-1 

0 

±1 

0 

0 


-1 

0 

0 

±1 

0 


-1 

0 

0 

0 

±1 

2 

+1 

±1 

0 

0 

0 


+1 

0 

±1 

0 

0 


+1 

0 

0 

±1 

0 


+1 

0 

0 

0 

±1 

3 

0 

0 

-1 

-1 

0 


0 

0 

1 

-1 

0 


0 

0 

-1 

1 

0 


0 

0 

1 

1 

0 


0 

-1 

0 

0 

-1 


0 

1 

0 

0 

-1 


0 

-1 

0 

0 

1 


0 

1 

0 

0 

1 


0 

0 

0 

-1 

-1 


0 

0 

0 

1 

-1 


0 

0 

0 

-1 

1 


0 

0 

0 

1 

1 


0 

-1 

-1 

0 

0 


0 

1 

-1 

0 

0 


0 

-1 

1 

0 

0 


0 

1 

1 

0 

0 


0 

0 

-1 

0 

-1 


0 

0 

1 

0 

-1 


0 

0 

-1 

0 

1 


0 

0 

1 

0 

1 


0 

-1 

0 

-1 

0 


0 

1 

0 

-1 

0 


0 

-1 

0 

1 

0 


0 

1 

0 

1 

0 

4 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 
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disturbing the equivalence property of the design. 


To define K for an MWP Method BBD, we set V$ = 0 in Equation (3). Based on this form 
of K and Equation (4), the conditions that must be satisfied, beyond those related to the 
HTC factors that are satisfied by the split-plot structure, are 

In^n.SQ,. = l^m'o + WcyMw 

rql n = n 0 l ni + W Q .u w V i. (10) 

The first condition is unchanged from the balanced case where m 0 and are defined 
in Table 1. The second is similar to the VKM unbalanced case, except that there are no 
conditions on the subplot pure quadratic terms. Recall the form of uw is Al w , and with 
w = 1 we have u w as a scalar equal to A. We rewrite the second condition in Equation (10) 
for w = 1 as 


in - • AWq,, V i, j , 

where Wq is a n x 1 vector for w = 1 of the form z\ l ni . Therefore, for a BBD we only 
need to consider the two possible values of z\ as 0 and 1. Table 9 summarizes the two cases 
of whole-plot types. For z\ = 0, we see that n 0 = rq = 2 k(k — 1) + n c . For z\ = 1, we End 

A = n 2 - n 0 
A = 2 k — [2 k(k — 1) + n c ] 

A = — ( 2k 2 — 4 k + n c ) . 

Consider a design for D(l,3) in Table 10. This design was derived from a balanced design 
with m — 3, n — 13, and N = 39. The unbalanced version maintains 3 whole-plots, however 
the overall size has been reduced to N — 26, by removing the replicated design points in the 
first two whole-plots. This is a similar construction as the design in Table 8 except we have 
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Table 9: MWP Method BBD Whole-Plot Types. 


Case 

w q, 7 

n t 

1 

0 

n i = 2 k{k — 1) + n c 

2 

1 

n 2 = 2 k 


now been able to include the subplot centers in whole-plot number 3, instead of requiring 
an additional whole-plot. For the design in Table 10, K is defined by no = 14, % = 6, 
u w = —8, rn'o = 81*,, M w = —614, and V s = 0. 

Note that a completely randomized version of a BBD with / = 4 contains N = 24 + n c 
design points. Therefore, the equivalent estimation design in Table 10 achieves the minimum 
number of whole-plots to support a second-order model and requires no additional subplot 
runs over the completely randomized version. In fact, this is a general result for MWP 
method BBD’s for D(l,/c), k G {2,3,4}. 

A summary table of unbalanced MWP method BBD’s is provided in Table 5. To achieve 
the minimum number of whole-plots, these designs require a large number of subplot runs 
per whole-plot, similar to the VKM method BBD’s. 

Minimal Point Designs 

As previously discussed, in some experimental situations the cost of the whole-plots ap- 
proaches that of the subplots and both are expensive. In this situation, it is desirable to 
simultaneously minimize the number of whole-plots and the total number of subplot runs. 
Therefore, a saturated or near-saturated designs might be considered. As in the CRD con- 
text, we recommend the use of the CCD or BBD, unless they are cost prohibitive. In 
addition, we recommend a cost comparison to conducting a completely randomized experi- 
ment, as proposed by Bisgaard (2000), which provides increased precision of the whole-plot 
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Table 10: MWP Method Unbalanced Box-Behnken Design for D(l,3). 
(z\ is a HTC factor, xi, X 2 , and X 3 are ETC factors.) 


Whole-Plot Zi 

Xi 


X3 

1 

-1 

-1 

0 

0 


-1 

1 

0 

0 


-1 

0 

-1 

0 


-1 

0 

1 

0 


-1 

0 

0 

-1 


-1 

0 

0 

1 

2 

1 

-1 

0 

0 


1 

1 

0 

0 


1 

0 

-1 

0 


1 

0 

1 

0 


1 

0 

0 

-1 


1 

0 

0 

1 

3 

0 

0 

-1 

-1 


0 

0 

-1 

1 


0 

0 

1 

-1 


0 

0 

1 

1 


0 

-1 

-1 

0 


0 

-1 

1 

0 


0 

1 

-1 

0 


0 

1 

1 

0 


0 

-1 

0 

-1 


0 

-1 

0 

1 


0 

1 

0 

-1 


0 

1 

0 

1 


0 

0 

0 

0 


0 

0 

0 

0 
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factors, versus the split-plot structure. 


Goos (2002) highlights that all saturated designs are equivalent. For a saturated design, we 
have N design points to estimate a model with p parameters where N = p. Therefore, the 
model matrix, X is a square full rank matrix. The OLS estimate is 

Pols = (X'X)- 1 X'y = x'-^X'r'X'y = X^y. 

The GLS estimate is 

Pols = (X'E-'XJ^X'E-V = X- 1 S(X')" 1 X / S~ 1 y = X^y, 

which is the same as the OLS estimate. Even though all saturated designs are equivalent, 
some are more suitable for augmentation to provide pnre-error estimates of the variance 
components and lack-of-fit degrees of freedom while maintaining the equivalence property. 

In the previous sections, design construction was based on an analytical derivation of the 
form of the K matrix to satisfy the equivalence condition for a particular design structure. 
As previously mentioned, Equation (2) can also be used to numerically test the equivalence 
of any arbitrary split-plot design without knowledge of the form of K. If the design is 
found to be equivalent, then K can be computed directly from Equation (1). Computing 
K provides insight into permissible augmentation strategies that maintain the equivalence 
property. For example, we can determine if the design can be augmented by adding subplot 
centers one-by-one or a whole-plot containing all factors at their center level. 

Many of these minimal point designs feature non-orthogonal subplot designs. While orthog- 
onality is desirable and was maintained in the CCD and BBD constructions, these designs 
demonstrate the flexibility in generating equivalent estimation designs. In addition, to ob- 
tain minimal point split-plot designs, we have relaxed our constraint of n mm = 2, allowing 
whole-plots of size one, which are typically undesirable for HTC factors. 
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Tabic 11: Notz Saturated Design Design for D(l,3). 
( Zi is a HTC factor, aq, x 2 , and x 3 are ETC factors.) 


Whole-Plot 

Zl 

X \ 

X 2 


1 

-1 

1 

1 

-1 


-1 

1 

-1 

1 


-1 

-1 

1 

1 


-1 

1 

1 

1 


-1 

-1 

-1 

-1 

2 

1 

1 

-1 

-1 


1 

-1 

1 

-1 


1 

-1 

-1 

1 


1 

1 

1 

-1 


1 

1 

-1 

1 


1 

-1 

1 

1 


1 

0 

0 

0 

3 

0 

1 

0 

0 


0 

0 

1 

0 


0 

0 

0 

1 


In this section, we consider the minimal point designs proposed by Notz (1982) to illustrate 
the utility of the numerical verification procedure. Consider a saturated Notz design for 
D(l,3) in Table 11. This design can be augmented to provide 1 degree of freedom for 
laek-of-fit ( dfLoF ), by including a subplot center run to whole-plot number 3, or by adding 
a subplot run to whole-plot number 2 with all of the subplot factors at their low level, 
[ — 1 —1 — 1 ]. Both of these augmentations can be verified to possess the equivalence 
property using Equation (2). A similar procedure can be employed for other numbers of 
factors. When the augmentation is not intuitive, a D-optimal augmentation strategy was 
applied to add the minimum number of design points while maintaining the equivalence 
property. Table 12 contains a summary of Notz equivalent estimation split-plot designs. 

Following a similar approach, equivalent estimation split-plot versions of the Hoke (1974), 
Box and Draper (1974), hybrid proposed by Roquemore (1976), and small composite pro- 
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Table 12: Summary of Notz Equivalent Estimation Split-Plot Designs. 

(w is the number of HTC factors, k is the number of ETC factors, m is the number of 
whole-plots, n is the number of subplot runs per whole-plot denoted by the whole-plot size 
with its frequency in parentheses, N is the total number of subplot runs, and c//l 0 f is the 

number of lack-of-fit degrees of freedom.) 


Factors Saturated Augmented 


w 

k 

m 

n (f'req. ) 

N 

m 

n( f'req.) 

N 

df LoF 

1 

1 

3 

3(1), 2(1), 1(1) 

6 

3 

3(1), 2(2) 

7 

1 

1 

2 

3 

5(1), 3(1), 2(1) 

10 

3 

5(1), 4(1), 2(1) 

11 

1 

1 

3 

3 

7(1), 5(1), 3(1) 

15 

3 

7(1), 5(1), 4(1) 

16 

1 

1 

4 

3 

11(1), 8(1), 2(1) 

21 





2 

1 

7 

2(3), 1(4) 

10 

7 

2(4), 1(3) 

11 

1 

2 

2 

7 

3(3), 2(2), 1(2) 

15 

7 

4(1), 3(2), 2(2), 1(2) 

16 

1 

2 

3 

7 

6(1), 4(3), 1(3) 

21 





2 

4 

7 

9(1), 5(3), 2(1), 1(2) 

28 

7 

9(1), 5(3), 3(1), 2(1), 1(1) 

30 

2 

3 

1 

12 

2(3), 1(9) 

15 

12 

2(8), 1(4) 

20 

5 

3 

2 

12 

2(9), 1(3) 

21 





3 

3 

12 

4(1), 3(4), 2(5), 1(2) 

28 
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posed by Hartley (1959), Draper (1985) and Draper and Lin (1990) minimal point response 
surface designs have been constructed. 

The ability to numerically verify that a design provides equivalent estimates using Equa- 
tion (2) has proven useful in exploring other families of RSM designs that do not provide 
a consistent analytical form of K, for example classes of D-optimal designs. In addition, 
numerical verification provides a convenient method to determine if model projection has 
disturbed the equivalence property. 

Concluding Remarks 

In this paper, we have relaxed the assumption of balance and demonstrated that unbalanced 
equivalent estimation split-plot designs can be constructed from families of completely ran- 
domized response surface designs. We have shown how to start with a balanced design and 
reduce the size of the whole-plots that contain subplot replicates. These unbalanced designs 
are applicable when the subplot cost becomes relatively expensive, approaching the cost of 
the whole-plots, and the experimental apparatus is suitable for unequal whole-plot sizes. 

The construction methods discussed in this paper provide a flexible approach to transform 
completely randomized response surface designs into a split-plot structure that provides best 
linear unbiased estimates and simplified estimation. A comprehensive catalog containing 
split-plot versions of many classical response surface designs is available at the author’s 
website, and includes balanced and unbalanced versions. The straightforward and intuitive 
nature of these construction methods coupled with the pre-built catalog makes equivalent 
estimation split-plot designs an attractive option for practitioners. 
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